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Abstract 

A binary lattice gas model that allows for multiple occupancy of lattice sites, inspired by recent 
coarse-grained descriptions of solutions of interacting polymers, is investigated by combining the 
steepest descent approximation with an exploration of the multidimensional energy landscape, and 
by Gibbs ensemble Monte Carlo simulations. The one-component version of the model, involving 
on site and nearest neighbour interactions, is shown to exhibit microphase separation into two 
sub-lattices with different mean occupation numbers. The symmetric two-component version of 
the multiple occupancy lattice gas is shown to exhibit a demixing transition into two phases above 
a critical mean occupation number. 
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I. INTRODUCTION 



Simple fluids are dominated by excluded volume effects, which, within lattice gas models, 
are accounted for by the single occupancy constraint, whereby each site on a lattice can 
be occupied by at most one molecule. Effective interactions between macromolecules or 
self-assembled aggregates in complex fluids, on the other hand, can be "soft", i.e. lack an 
impenetrable core. A good example is the effective pair potential between the centres of 
mass (CM) of interacting polymer coils, obtained by taking statistical averages over monomer 
conformations for fixed distances r between the CMs l[ m ■ Recent extensive simulations of 
self-avoiding walk polymers carried out over a wide range of concentrations show that the 
repulsive state-dependent effective CM pair potential is of roughly Gaussian shape, of width 
overned by the polymer radius of gyration, and of maximum amplitude v (r = 0) ~ fc^T 
3|; this behaviour reflects the fractal nature of polymers in good solvent, leading to a 
low entropic cost at full overlap. Several other complex fluids have been shown to exhibit 
nltrasoft coarse-grained interactions,^, star po.ymers Q or effective particles considered 
within dissipative particle dynamics p. 

The penetrant of the correspond "Ganssran core" (GC) mode,, ^l^(-r% 
leads to interesting phase behaviour at low temperatures (T* = ksT/e <C 1), [a,l2|, but under 
conditions relevant for polymer solutions (T* « 1), the model behaves like a "mean field" 
fluid 81. However binary Gaussian core mixtures, characterised by different energy scales e n R 

nun 

for the three types of pairs, lead to phase separation for moderate couplings (8J, 13, 110] • This 
demixing, which occurs for purely repulsive, penetrable interactions, is of a very different 
nature than the usual phase separation of incompatible fluids, which is generally driven by 
the long-range attractive intermolecular forces. 

In this paper we examine the simplest lattice gas representation of penetrable particles, 
in an effort to gain further insight into this novel class of phase transitions. The penetrable 
nature of the effective interaction is reflected by allowing multiple occupancy of each lattice 
site. The simplest version of the model involves only on-site interactions of particles of the 
two species, with different energy penalties for the different types of pairs. While a steepest 
descent estimate of the grand partition function predicts phase separation, it will be shown 
that a more accurate treatment reveals the spurious nature of this transition, as expected for 
an effectively zero- dimensional system. The addition of nearest neighbour interactions leads 



2 



to a genuine demixing transition, similar to that predicted for continuous versions of soft 
core mixtures examined earlier within the random phase approximation (RPA) j^, . We 
explore this behaviour by mean field theories and Gibbs ensemble Monte Carlo simulations. 
We also study the topology of the energy landscape. 



II. THE MULTI-OCCUPANCY MODEL 



Consider a (^-dimensional lattice of L sites and coordination number q, which can ac- 
commodate particles of two species. The occupancy of each site is characterised by the 
two-component vector of integer occupation numbers rij = (n[ , 1 ), 1 < i < L. Particles 
interact only when they are on the same site or on nearest neighbour (n.n.) sites. The 
on-site and off-site couplings are characterised by 2 x 2 matrices of interaction energies e a p 
and r] af 3 (1 < a, (3 < 2). In terms of the occupation numbers, the energy of the system reads: 
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where the last summation is over n.n. sites only. For given chemical potentials fi = (//i,//2 
of the two species, the grand partition function is: 
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The first sum is over all possible distributions of occupation numbers of the L sites, and the 
effective chemical potentials //* = // a — e aa /2, /j,* = (/i*,/^) have been introduced. 

The main task of this paper is to approximately evaluate the grand partition function (J2J) 
as a function of the temperature T {(3 = l/{ksT)) and of the chemical potentials [x\ and \x\ 
of the two species. Using the standard rules of Statistical Mechanics, this will lead directly 
to the pressure P = /csTln(H)/U and to the mean occupation numbers n\ and n\, and hence 
to the composition of the mixture. Phase separation will be signalled by a discontinuous 
change of the mean occupation number for well-defined values of P, and \x\. We shall 
first consider the simplest version of the model which only involves on-site interactions (i.e. 
r) = 0), before examining the case including nearest neighbour interactions. 



III. THE CASE OF ON-SITE INTERACTIONS ONLY 



In the absence of couplings between particles on different sites, the grand partition func- 
tion (J2J), with 77 = 0, factorises into L identical single site functions S site : 
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The single site partition sum is not tractable analytically. For sufficiently large chemical 
potentials, implying mean occupation numbers n^nlj ^> 1, the latter may be treated as 
continuous variables, and the sums in eq. (j3J) replaced by integrals: 

exp {dyb* ■ n — in* ■ e ■ n) 
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The integrand is sharply peaked around the most probable occupation numbers n* = 
(77,1,712). Adopting the standard steepest descent method, we expand the logarithm of the 
integrand around n* to second order in Sn = n — n*, i.e.: 
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where ip{n) = dlnr(n)/dn is the digamma function, and the covariance matrix a is given 
by: 



(6) 



with Tp^'in) = dip(n)/dn (trigamma function). The location of the maximum follows from 
the extremum condition 

PtS-P^-e- r K 1 J =0, (7) 

which provides the relation between /a* and n*. Inserting the expansion (jSJ) into the integrand 
in eq. (jlj), and extending the integral to negative occupation numbers (which is again justified 



provided n* 3> 1), the resulting Gaussian integral is easily evaluated and leads directly to 
the following equation of state: 

(3Pv = n* + n* + In*' ■ e ■ n* + In I -jL= ] , (8) 
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where vq = V/L is the volume per site. 

The virial pressure in eq. (|5|). derived on the basis of steepest descent, is thermodynam- 
ically inconsistent with the relation (|7J) between chemical potentials and mean occupation 
numbers (identical to their most probable values). If the latter is integrated, an expression 
for the pressure follows which coincides with eq. (jHJ), but without the logarithmic fluctua- 
tion term. The situation is reminiscent of the RPA treatment of continuous penetrable (e.g. 
Gaussian) core mixtures, where the virial and compressibility routes lead to equations of 
state similar in structure to eq. (jHJ) with (virial) and without (compressibility) the last term 
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Il0l | . The two equations of state become identical in the high density limit pointing 
to the asymptotic "mean-field fluid" nature of systems of penetrable particles. 

In view of the above analogy, it is not surprising that the equation of state (jEJ) 
leads to phase separation between sites with different mean occupation numbers when 
e i2 > e n e 22- This predicted phase separation is obviously spurious, because of the effec- 
tively O-dimensional nature of the model when restriction is made to on-site couplings. The 
above steepest descent treatment has two obvious shortcomings. The first deficiency is the 
replacement of the discrete sums in eq. Q by an integral, as well as the extension of the 
lower bounds in the integral to — oo; these approximations are expected to break down 
for low mean occupation numbers n* , or equivalently for low chemical potentials /z* . This 
deficiency can be easily removed by reverting from the integrals to discrete sums over an 
interval of values of n* in the vicinity of the most probable values n*, and retaining the con- 
tinuous integrals outside this interval. In the one-component case, the resulting single-site 
partition function reads: 
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where ri\ ow = max( \n*~\ — m, 0) and nhigh — 1 = + m — 1, and the ceil function |~] rounds 
to the nearest higher integer. 



The steepest descent result (JHJ) is recovered for m = and sufficiently high n*, such that 
erf(n*A/cr/2) = 1. Results for the equation of state are shown in Fig. [T] Eq. (jHJ) is seen to 
yield surprisingly accurate results, even at relatively low densities (n* ~ 1), and only breaks 
down in the n* — > limit. 




FIG. 1: Equation of state of the one-component lattice gas with on site interaction energy e = 
lksT. The normalised pressure (3Pvq is plotted as a function of the mean occupation number 
(n) of one cell. Shown are values obtained by using the steepest descent method (dashed curve, 
circles), the pressure consistent with the chemical potential (J7J) (dash-dotted line, squares), and 
the pressure obtained from the grand canonical partition sum by explicitly summing 2m terms 
around the maximum of the Boltzmann factor as explained in the text with m = 100 (continuous 
line, triangles). At this value the correction term in equation Q is well below 10 -5 , so that the 
resulting pressure can be considered nearly exact. The steepest descent method is seen to be very 
accurate except, as expected, at low densities. 

A more fundamental shortcoming of the steepest descent method is that eq. Q in fact 
predicts two maxima n*, separated by a saddle point. The situation is pictured in Fig. |21 
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The steepest descent method includes only the contribution of the highest (global) maximum 
of the integrand to the partition function (jlj). This approximation leads to a discontinuous 
jump of the composition vector n* when the two maxima become equal, and hence to a 
first order phase transition. In reality, when the contributions of the two maxima of the 
configuration integral (J3J) are properly included, the mean occupation numbers (ni) and 
(712) vary continuously with the chemical potentials of the two species, as illustrated in 
Fig. El according to expectation. 




FIG. 2: Boltzmann factor exp(/3/^ • n — f3E{rv)) / '(niln^!) as a function of the mean occupation 
numbers n\,n2 for the lattice system with on site interaction c\\ — 622 — 612 — 

and no nearest neighbour interaction, at the chemical potentials /ii = 6.5/cbT, [12 = 6.0ksT. The 
appearance of a secondary maximum is clearly seen. 
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FIG. 3: Mean occupation numbers (wi), (^2} in a symmetric binary lattice gas mixture for on site 
interaction energies en = €22 = 0.1&bT €12 = 0.2/cbT and vanishing nearest neighbour interaction 
as a function of chemical potential fi2, keeping /j>i = 6.5/c#T fixed. The values were calculated using 
a simple generalization of equation Q As expected, the occupation numbers vary continuously with 
chemical potentials in this effectively 0-dimensional system. 

IV. INCLUDING NEAREST NEIGHBOUR INTERACTIONS 

We now consider the model denned by eqns. and (0) and including n.n. interactions 
between particles of both species, i.e. 77 7^ 0. The partition function (j2J) no longer factorises 
into single site sums. Proceeding as in Section ITTT1 we replace the discrete sums over occu- 
pation numbers by integrals and approximate the integrand by Gaussians centred on each 
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of the local maxima. The latter are determined by a generalisation of eq. (|7|). namely 



(3^ = (3nf-e+J2^T- r J + 



1>(<i + 1) 



1 < i < L. 



(10) 



In the limit 77 — > 0, the n* on different sites are independent and can each take two val- 
ues, corresponding to the two maxima found in Section 11111 When the n.n. interactions 
are switched on, the n* are coupled so that a complex energy landscape emerges in 2L- 
dimensional occupation-number space, featuring a rapidly increasing number of maxima of 
the Boltzmann factor. By continuity one can still expect 2 L solutions to eqns. (fTUj) in the 
most general case, at least for sufficiently weak coupling. The approximate grand partition 
function now reads 



S = Yl exp 



(ij) 
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where the covariance matrices <Tj are obvious generalisations of eq. Q, with n* replaced by 



The grand partition function is thus expressed as a sum over all "configurations" {n*} 
which maximise the integrand in the continuous representation generalising eq. (j3J), i.e. which 
satisfy the set of eqns. (|1U|). The Gaussian fluctuation integrals in the second factor depend 
on the configuration {n*} through the covariance matrices ov For given {er^} the coupled 
Gaussian integrals can be calculated explicitly by a standard normal mode analysis. 

The locations of the local maxima n* depend on the occupation numbers of the neigh- 
bouring cells n* via (|10p. Were that not the case, then the partition function (II 1J) would 
exactly factor into a product of a partition sum that is isomorphic to that of an Ising model 
(first factor in (fTTj) ) and the partition sum of a Gaussian model. The phase behaviour of 
a very similar system has been investigated in Ref. 11]. The model we examine can thus 
only be mapped onto an Ising (plus Gaussian) model in the limit of weak nearest neighbour 
interactions rj. 

In order to gain further insight, we consider the "ground state", corresponding to the 
low temperature or strong coupling limit. In that limit all n* corresponding to the lowest 



minimum on the total energy surface are expected to be equal (n* = n*), and determined 
by a single extremum condition: 



PfS = pn* 1 ■ (e + qrj) + 



v(nj + 1; 
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where q is the coordination number of the lattice. In the low temperature limit only the 
contribution of this one configuration to the sum in equation (|11|) is taken into account. The 
second factor in this equation then becomes invariant with respect to translations, so that 
the corresponding integral is best evaluated in Fourier space. 

We therefore introduce the Fourier components of the density fluctuations 



(13) 



where denotes the position of the lattice site j in units of the lattice constant. Assuming 
Born-von Karman boundary conditions the wave vector k can only assume discrete values. 
The total interaction energy can be written in Fourier space as 



(14) 
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where the matrix D(k) is given by 
D(k) = o- + /3r7 eik ' Rj - 

j'n.n.O 

For a simple cubic lattice we obtain 



(15) 



D(k) = <x + 2 (cos k x + cos k y + cos k z )/3rj. 



(16) 



The Fourier modes of the density fluctuations thus decouple, so that the multiple integral 
in equation (JTTJ) can be evaluated in Fourier space as a product of Gaussian integrals 



dSnx ■ ■ ■ ddriL exp 



L I ddh^ • . . / d5n^ L exp 



(k)| 
(17) 
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The grand partition function then reduces to: 

exp | in* • (e + qr)) ■ n* + £}<V« a + 1) - lnT« + 1)1 

x(27r) i exp|-||^ln(|D(k)|)} ! (18) 

leading to the equation of state 

^ = ^KV«+l)-lnr« + l)] + in*.(e+gr7)-n*+l n (27r)-i J -^ln (|D(k)|) . 

(19) 

In the limit of large occupation numbers (n* 3> 1) and vanishing nearest neighbour interac- 
tion (rj = 0) we recover equation (JHJ). As in the case without nearest neighbour interaction, 
there is a thermodynamic inconsistency between the e.o.s. (fT9*|) and the dependence of the 
chemical potentials on the occupation numbers (|T2|) . Integrating the latter would again lead 
to an e.o.s. similar in structure to equation ()19|). albeit without the final two logarithmic 
terms. The virial equation of state is shown in Figure EJ For high nearest neighbour interac- 
tion (e — qrj no longer being positive definite) the fluctuation term diverges at k = (n, 7T, n). 
This instability indicates a transition of the ground state from the homogeneous phase to a 
microfluid phase where high and low occupation numbers alternate in a checkered fashion. 
To show that the microphase exists also for the one-component fluid (for qr] > e) a canonical 
Monte Carlo simulation was run on a 10 x 10 x 10 cubic lattice (q = 6) with the parameters 
e = 0.1/cbT, t] = 0.033ksT. Initially every lattice site was filled with 20 particles. After 
equilibration a histogram of the occupation numbers occuring over the next 1 000 000 steps 
was recorded. The result is shown in Figure It is clearly seen that the initial peak at 
n = 20 has vanished and a bimodal distribiution around n pa and n pa 40 has developed. 
The development of the microphase is easy to understand when we estimate and compare 
the energetic and entropic costs of the microphase with (mean) occupation numbers n a and 
rib with those of a homogeneous phase with average density n = (n a + n^)/2. The entropy of 
the microphase per site is s micro = ks[n a lnri a + nblnnf,]/2, compared to that of the homoge- 
neous phase s hom - = A; B nlnn. The difference between the two is simply the mixing entropy 
of the ideal gas, thereby favouring the homogeneous density distribution. On the other hand 
the internal energy of the microphase per site is ^ micro = + + qr]n a nf ) /2, com- 

pared to that of the homogeneous phase M hom - = (e + qr])n 2 /2. The difference between the 
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two internal energies is ■u micro — -u hom = (e — qrf)(n b — n a ) 2 /8. For strong nearest neighbour 
interaction the microphase is therefore energetically favoured. The nature of the transition 
from the homogeneous state to the micro state will be investigated elsewhere. 




n 

FIG. 4: Virial equation of state of a one component lattice system with on site interaction e = lfceT 
and nearest neighbour interaction rj = O.lksT (triangles), 77 = 0.2fc^T (circles), and rj = 0.3/cbT 
(squares). For (3rj > 1/6 our analysis predicts an instability of the homogeneous phase towards 
microphase separation beyond a certain density. For rj = O.lksT this criterion is not obeyed and 
the fluid remains stable at all densities. For rj = 0.2/cbT the instability is at n* = 4.48 (outside the 
density range of this plot), and for rj = 0.3 it is at n* = 0.69. 

The explicit knowledge of the Fourier components of the density fluctuations allows us 
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FIG. 5: Histogram of the occupation numbers in a one component lattice system (10 x 10 x 10) with 
e = O.lksT and r/ = 0.033A;£T, The probability p(n) to find a lattice site with occupation number 
n (averaged over 1 000 000 steps after equilibration) is shown. The canonical simulation was 
started with 20 particles on each lattice site. Because of the large nearest neighbour interaction 
the initially homogeneous fluid develops into a microphase with a bimodal distribution of the 
occupation numbers. 

to evaluate the structure factor of the fluid. For we obtain 

S a p{k) = (8n a (k)5np{-k)) 

1 / dSrn(k) J d5n 2 {k)5n a {k)S~n {k) exp(-<fn(k) • D(k) ■ <fn(k)) 
" n\+n* 2 ] d5~m{k) j d5n 2 (k) exp(-<5~n(k) ■ D(k) ■ <5n(k)) 

From the particle-particle structure factors S a p the Bhatia-Thornton structure factors 
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can be extracted. These are compared with data from a canonical Monte Carlo simulation 
in Figure Elfor the phase point (ni) = 10, (712) = 1. The agreement is reasonable, although 
there are still fairly large statistical errors, a problem commonly encountered in direct calcu- 
lations of S(k) (as opposed to Fourier transforming g(r)). The error bars on the simulation 
results (not shown) are rather large. 




FIG. 6: Bhatia-Thornton structure factors for the symmetric binary mixture with (3e\\ = $€22 = 
0.1, (3e\2 = 0.2, ij = O.le. The mean-field values (open symbols, lines) are compared to simulation 
data of a canonical Monte Carlo simulation (filled symbols) on a 10 x 10 x 10 lattice. Shown 
are the density-density (S nn , triangles), density-concentration (S nc , circles) and concentration- 
concentration (S cc , squares) structure factors along the Ill-direction for the occupation numbers 
(ni) = 10, (712) = 1. The statistics for the simulation data is known to be poor in A;-space. 

In the calculations above we neglected all contributions from configurations in the par- 
tition sum (jllj) other than the (presumed) ground state, where all occupation numbers are 
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the same. This is of course only justified when the other configurations have much higher 
energies. To check this hypothesis, we explore the energy landscape of small lattice systems 
with fii = /i 2 . For a 2 x 2 x 2 lattice systems, one can explicitly enumerate all local energy 
minima. In the case of no nearest neighbour interactions the Boltzmann factor has 2 8 degen- 
erate local maxima. Slowly switching on nearest neighbour interactions by setting 77 = Ae 
and varying A, one can numerically follow the local maxima of the Boltzmann factor. Due 
to the high symmetry, many local minima stay degenerate even at finite nearest neighbour 
interaction. When increasing 77, some local minima turn into saddle points and vanish, as 
illustrated in Figure [H The remaining excited energy states depending on A are shown in 
the same figure. The increasing gap between the ground state and the excited states is 
clearly seen. For larger lattices the explicit enumeration of all local minima is no longer 
feasible because of their large number. We therefore resort to a simple Metropolis algorithm 
employed in numerical optimization problems Q|. Instead of systematically searching the 
2L dimensional occupation number state, we generate trial configurations randomly from 
the previously found local minimum, analogous to the random particle displacements used 
in Monte Carlo simulations. After each random step in the space of occupation numbers a 
local minimization of the energy is performed. The random step is then accepted or rejected 
based on the energy difference of the two local minima. Effectively the very rugged energy 
landscape is therefore replaced with a locally constant function, thereby eliminating the 
numerical problems one usually encounters when solving global optimization problems. In 
Figure |H1 the resulting local energy minima are shown together with the (presumed) homo- 
geneous ground state. All the minima lie on or slightly above the presumed ground state, 
thereby supporting our hypothesis that the ground state is the homogeneous state. In Fig- 
ure |H1 the difference between the local minima and the ground state energy is plotted. Even 
though some local minima are not found during the search, one can clearly see a trend of 
an increasing gap between the ground state and the first excited state when the nearest 
neighbour interaction increases. For a nearest neighour interaction that equals 10 percent 
of the on-site interaction one finds a gap of approximately 1.5/cbT. Even stronger nearest 
neighbour interaction should further widen the gap, leading to a higher accuracy of the 
ground state approximation. 
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FIG. 7: Local energy minima for a 2 x 2 x 2 lattice system at [i\ = lii = 6.5/cbT. The energy 
minima are degenerate for vanishing nearest neighbour interaction (A = 0) and split for finite A 
into sets of local minima related by symmetry of the lattice. At certain critical A, some minima 
turn into saddle points and vanish, as visualized in the inset. 
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FIG. 8: Absolute ground state energies in a 5 x 5 x 5 lattice at /3fj,\ = (3^2 = 6.5 as found by the 
optimization procedure described in the text (plusses) compared to the energy of the homogeneous 
state, as a function of the nearest neighbour interaction r] = Xe. On site interactions are en = 
e 22 = O.lfcfiT, ei2 = 0.2k B T. 
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FIG. 9: Local energy minima with respect to the ground state, as found by the optimization 
procedure described in the text. The on site interactions are en = £22 = O.lksT, eu = 0.2fceT. 
The nearest neighbour interactions are rj = Ae with A = 0,0.005,0.01, . . . ,0.1. It is clearly seen 
that the energy gap from the ground state to the first excited state increases with increasing nearest 
neighbour interaction. The ground state approximation should therefore become more accurate 
for stronger nearest neighbour interaction. 
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V. GIBBS ENSEMBLE SIMULATIONS 



To test the validity of the ground state approximation, Gibbs ensemble Monte Carlo 
simulations w| of the model have been performed on a 10 x 10 x 10 lattice. After an 
equilibration time of 1 000 000 steps two lattices were allowed to exchange particles and 
random moves within each lattice were attempted according to the standard Metropolis 
algorithm. The simulations were restricted to symmetric interactions with e n = €22 and 
^11 = f]22- Since in these systems the equation of state must be invariant with respect to 
relabelling of the species (or equivalently mirroring the concentration x — > 1 — x), volume 
changes of the lattice were not necessary, i. e. the two lattice systems were automatically 
under equal pressure. Two typical phase diagrams are shown in Figure ^] and ^2 It is seen 
that the ground state approximation predicts the correct shape of the phase diagram, albeit 
underestimates the densities of the coexistencing phases. The discrepancy becomes larger 
for larger interaction energies, when the contributions of the excited states are no longer 
negligible, and the ground state approximation becomes less accurate. 
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FIG. 10: Phase diagram of a symmetric binary lattice system with en = £22 = 0.1&bT, €12 = 
0.2/cbT, r) = O.le. Shown are the phase coexistence points of the Gibbs ensemble simulation 
(crosses) and the predictions of the binodal (solid line) and spinodal (dashed line) based on the 
ground state approximation. No microphase separation is expected for such low nearest neighbour 
interactions. 
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. 11: Phase diagram of a symmetric binary lattice system with en = 622 = 0.5/cbT, €12 = lksT, 
O.le as in Figure ITU1 
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VI. CONCLUSION 



The lattice model introduced in this paper, which allows for multiple occupancy of lattice 
sites by one or two species, with positive energy penalties, and includes repulsive nearest 
neighbour interactions, leads to intersting phase behaviour. In the one-component version, a 
steepest descent analysis allowing for gaussian fluctuations, predicts microphase separation 
onto two interpenetrating sub- lattices with different mean occupation numbers (" checkered" 
phase) when the ratio r]/e of nearest-neighbour to on-site couplings exceeds a threshold. This 
prediction is confirmed by Monte Carlo simulations. 

The two-component extension of the model leads to a macroscopic demixing transition 
into phases of different concentrations of the two species. The steepest descent analysis 
predicts such a phase separation, even in the absence of nearest-neighbour interactions, 
but this is shown to be due to the expected break-down of the approximation. In the 
presence of nearest-neighbour interactions, the steepest-descent analysis corresponds to a 
"ground state" approximation. A numerical search of local minima of the energy surface 
points to a significant gap (> lk B T) between the homogeneous "ground-state" and the 
first "excited states", thus providing support for the validity of the approximation. The 
predicted phase diagrams agree reasonably well with the results of Gibbs ensemble MC 
simulations. The present model is directly inspired by the continuous "Gaussian core" 
model, first introduced by Stillinger and which has recently been shown to provide an 
adequate coarse-grained description of interacting, interpenetrating polymer coils 3]. The 
corresponding lattice model, considered here, appears to have some unorthodox features. 
Contrary to more familiar lattice gas models, it is not isomorphous to an Ising spin model. 
Hence it is not yet clear whether its critical behaviour belongs to the Ising universality class. 
The binary version is expected to exhibit both the microphase separation found here in 
the one-component case, and the macroscopic phase separation anologuous to that of the 
continuous binary gaussian core model P, 0, The interplay between these two phase 
transitions should lead to interesting and novel behaviour, which will be explored in future 
work. 
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